################################################
## 17_combining_chains_all_processes.R
################################################
## Use this file to combine estimated chains from
## Brazil 2MPE
## Brazil 1MPE
## Mexico 2MPE
## Mexico 1MPE
################################################

rm (list=ls())

print (version)

# platform       aarch64-apple-darwin20      
# arch           aarch64                     
# os             darwin20                    
# system         aarch64, darwin20           
# status                                     
# major          4                           
# minor          1.2                         
# year           2021                        
# month          11                          
# day            01                          
# svn rev        81115                       
# language       R                           
# version.string R version 4.1.2 (2021-11-01)
# nickname       Bird Hippie 


#### Combining scans from all samples (CJRpressure model)

library (runjags)
library (coda)
library (reshape2)
library (random)
library (scales)

#####################
#### Brazil 2MPE #### 
#####################
rm(list=ls())

# Set working directory
savePathRep <- "/all_data/"
setwd (savePathRep)

# Change names in the next two lines

all.processes <- c("brazil_competing_final_result_sample_1a.Rda",
                   "brazil_competing_final_result_sample_1b.Rda",
                   "brazil_competing_final_result_sample_2a.Rda",
                   "brazil_competing_final_result_sample_2b.Rda",
                   "brazil_competing_final_result_sample_3a.Rda",
                   "brazil_competing_final_result_sample_3b.Rda",
                   "brazil_competing_final_result_sample_4a.Rda",
                   "brazil_competing_final_result_sample_4b.Rda",
                   "brazil_competing_final_result_sample_5a.Rda",
                   "brazil_competing_final_result_sample_5b.Rda",
                   "brazil_competing_final_result_sample_6a.Rda",
                   "brazil_competing_final_result_sample_6b.Rda",
                   "brazil_competing_final_result_sample_7a.Rda",
                   "brazil_competing_final_result_sample_7b.Rda",
                   "brazil_competing_final_result_sample_8a.Rda",
                   "brazil_competing_final_result_sample_8b.Rda",
                   "brazil_competing_final_result_sample_9a.Rda",
                   "brazil_competing_final_result_sample_9b.Rda",
                   "brazil_competing_final_result_sample_10a.Rda",
                   "brazil_competing_final_result_sample_10b.Rda")
                   

max.gelman.r <- c()
allChains <- mcmc.list()
for (k in 1:10){
  chain <- list()
  which.chains <- grep (paste0("sample_",k,"[ab]"), all.processes)

  # print sample iteration
  print (k)
  
  for (i in 1:length(which.chains)){
    load (all.processes[which.chains[i]])
    chain[[i]] <- as.mcmc(chainsComp); rm (chainsComp)
    # print process iteration within sample
    print (i)
  }

  gelman.Comp <- gelman.diag (chain, multivariate=F)
  max.gelman <- max (gelman.Comp$psrf[,1], na.rm=T)
  print (max.gelman)
  max.gelman.r <- c(max.gelman.r, max.gelman)

  allChains[[k]] <- as.mcmc.list (chain)
}


# Add convergence check for 20 chains all together
# This step ascertains that different chains *based on different samples*
# converge more or less on the same parameter values
# Because the samples are different, gelman hat values are somewhat larger
Chains <- as.mcmc.list (list (
  allChains[[1]][[1]], allChains[[1]][[2]], 
  allChains[[2]][[1]], allChains[[2]][[2]], 
  allChains[[3]][[1]], allChains[[3]][[2]], 
  allChains[[4]][[1]], allChains[[4]][[2]], 
  allChains[[5]][[1]], allChains[[5]][[2]], 
  allChains[[6]][[1]], allChains[[6]][[2]], 
  allChains[[7]][[1]], allChains[[7]][[2]], 
  allChains[[8]][[1]], allChains[[8]][[2]], 
  allChains[[9]][[1]], allChains[[9]][[2]], 
  allChains[[10]][[1]], allChains[[10]][[2]]))

gelman.all <- gelman.diag (Chains, multivariate=F)
max.gelman <- max (gelman.all$psrf[,1], na.rm=T)
print (max.gelman)



theta <- delta <- beta <- alpha <- eta <- lambda <- c()
voteDis <- voteAgr <- voteUnc <- deviance <- c()
for (i in 1:length(allChains)){
  process.1 <- allChains[[i]][[1]]
  process.2 <- allChains[[i]][[2]]
  
  th1 <- process.1[,grep("theta", colnames(process.1))]
  dt1 <- process.1[,grep("delta", colnames(process.1))]
  bt1 <- process.1[,grep("beta", colnames(process.1))]
  et1 <- process.1[,grep("^eta", colnames(process.1))]
  al1 <- process.1[,grep("^alpha", colnames(process.1))]
  lm1 <- process.1[,grep("lambda", colnames(process.1))]
  dis1 <- process.1[,grep("VoteDisagree", colnames (process.1))]
  agr1 <- process.1[,grep("VoteAgree", colnames (process.1))]
  unc1 <- process.1[,grep("probVoteUncertain", colnames (process.1))]
  dev1 <- process.1[,grep("deviance", colnames (process.1))]

  th2 <- process.2[,grep("theta", colnames(process.2))]
  dt2 <- process.2[,grep("delta", colnames(process.2))]
  bt2 <- process.2[,grep("beta", colnames(process.2))]
  et2 <- process.2[,grep("^eta", colnames(process.2))]
  al2 <- process.2[,grep("^alpha", colnames(process.2))]
  lm2 <- process.2[,grep("lambda", colnames(process.2))]
  dis2 <- process.2[,grep("VoteDisagree", colnames (process.2))]
  agr2 <- process.2[,grep("VoteAgree", colnames (process.2))]
  unc2 <- process.2[,grep("probVoteUncertain", colnames (process.2))]
  dev2 <- process.2[,grep("deviance", colnames (process.2))]
  
  th <- rbind (th1, th2)
  dt <- rbind (dt1, dt2)
  bt <- rbind (bt1, bt2)
  et <- rbind (et1, et2)
  al <- rbind (al1, al2)
  lm <- rbind (lm1, lm2)
  disagree  <- c(dis1, dis2) 
  agree     <- c(agr1, agr2) 
  uncertain <- c(unc1, unc2) 
  dev  <- c(dev1, dev2)

  theta  <- rbind (theta, th)
  delta  <- rbind (delta, dt)
  beta   <- rbind (beta, bt)
  alpha  <- rbind (alpha, al)
  eta    <- rbind (eta, et)
  lambda <- rbind (lambda, lm)
  voteDis <- c(voteDis, disagree)
  voteAgr <- c(voteAgr, agree)
  voteUnc <- c(voteUnc, uncertain)
  deviance <- c(deviance, dev)
}

items2export <- list (theta, delta, alpha, beta, eta, lambda, voteDis, voteAgr, voteUnc, deviance)
names (items2export) <- c ("theta", "delta", "alpha", "beta", "eta", "lambda", "voteDis", "voteAgr", "voteUnc", "deviance")

# Uncomment next line to save file
# save (items2export, file="objects_2MPE_Brazil.Rdata")


#####################
#### Brazil 1MPE #### 
#####################
rm(list=ls())

# Set working directory
savePathRep <- "/all_data/"
setwd (savePathRep)

# Change names in the next two lines
all.processes <- c("brazil_pressure_final_result_sample_1a.Rda",
                   "brazil_pressure_final_result_sample_1b.Rda",
                   "brazil_pressure_final_result_sample_2a.Rda",
                   "brazil_pressure_final_result_sample_2b.Rda",
                   "brazil_pressure_final_result_sample_3a.Rda",
                   "brazil_pressure_final_result_sample_3b.Rda",
                   "brazil_pressure_final_result_sample_4a.Rda",
                   "brazil_pressure_final_result_sample_4b.Rda",
                   "brazil_pressure_final_result_sample_5a.Rda",
                   "brazil_pressure_final_result_sample_5b.Rda",
                   "brazil_pressure_final_result_sample_6a.Rda",
                   "brazil_pressure_final_result_sample_6b.Rda",
                   "brazil_pressure_final_result_sample_7a.Rda",
                   "brazil_pressure_final_result_sample_7b.Rda",
                   "brazil_pressure_final_result_sample_8a.Rda",
                   "brazil_pressure_final_result_sample_8b.Rda",
                   "brazil_pressure_final_result_sample_9a.Rda",
                   "brazil_pressure_final_result_sample_9b.Rda",
                   "brazil_pressure_final_result_sample_10a.Rda",
                   "brazil_pressure_final_result_sample_10b.Rda")

max.gelman.r <- c()
allChains <- mcmc.list()
for (k in 1:10){
  chain <- list()
  which.chains <- grep (paste0("sample_",k,"[ab]"), all.processes)
  
  # print sample iteration
  print (k)
  
  for (i in 1:length(which.chains)){
    load (all.processes[which.chains[i]])
    chain[[i]] <- as.mcmc(chainsPress); rm (chainsPress)
    # print process iteration within sample
    print (i)
  }
  
  gelman.Comp <- gelman.diag (chain, multivariate=F)
  max.gelman <- max (gelman.Comp$psrf[,1], na.rm=T)
  print (max.gelman)
  max.gelman.r <- c(max.gelman.r, max.gelman)
  
  allChains[[k]] <- as.mcmc.list (chain)
}


# Add convergence check for 20 chains all together
# This step ascertains that different chains *based on different samples*
# converge more or less on the same parameter values
Chains <- as.mcmc.list (list (
  allChains[[1]][[1]], allChains[[1]][[2]], 
  allChains[[2]][[1]], allChains[[2]][[2]], 
  allChains[[3]][[1]], allChains[[3]][[2]], 
  allChains[[4]][[1]], allChains[[4]][[2]], 
  allChains[[5]][[1]], allChains[[5]][[2]], 
  allChains[[6]][[1]], allChains[[6]][[2]], 
  allChains[[7]][[1]], allChains[[7]][[2]], 
  allChains[[8]][[1]], allChains[[8]][[2]], 
  allChains[[9]][[1]], allChains[[9]][[2]], 
  allChains[[10]][[1]], allChains[[10]][[2]]))

gelman.all <- gelman.diag (Chains, multivariate=F)
max.gelman <- max (gelman.all$psrf[,1], na.rm=T)
print (max.gelman)



theta <- delta <- beta <- alpha <- eta <- lambda <- deviance <- c()
for (i in 1:length(allChains)){
  process.1 <- allChains[[i]][[1]]
  process.2 <- allChains[[i]][[2]]
  
  th1 <- process.1[,grep("theta", colnames(process.1))]
  dt1 <- process.1[,grep("delta", colnames(process.1))]
  bt1 <- process.1[,grep("beta", colnames(process.1))]
  et1 <- process.1[,grep("^eta", colnames(process.1))]
  al1 <- process.1[,grep("^alpha", colnames(process.1))]
  lm1 <- process.1[,grep("lambda", colnames(process.1))]
  dev1 <- process.1[,grep("deviance", colnames (process.1))]
  
  th2 <- process.2[,grep("theta", colnames(process.2))]
  dt2 <- process.2[,grep("delta", colnames(process.2))]
  bt2 <- process.2[,grep("beta", colnames(process.2))]
  et2 <- process.2[,grep("^eta", colnames(process.2))]
  al2 <- process.2[,grep("^alpha", colnames(process.2))]
  lm2 <- process.2[,grep("lambda", colnames(process.2))]
  dev2 <- process.2[,grep("deviance", colnames (process.2))]
  
  th <- rbind (th1, th2)
  dt <- rbind (dt1, dt2)
  bt <- rbind (bt1, bt2)
  et <- rbind (et1, et2)
  al <- rbind (al1, al2)
  lm <- rbind (lm1, lm2)
  dev  <- c(dev1, dev2)
  
  theta  <- rbind (theta, th)
  delta  <- rbind (delta, dt)
  beta   <- rbind (beta, bt)
  alpha  <- rbind (alpha, al)
  eta    <- rbind (eta, et)
  lambda <- rbind (lambda, lm)
  deviance <- c(deviance, dev)
}

items2export <- list (theta, delta, alpha, beta, eta, lambda, deviance)
names (items2export) <- c ("theta", "delta", "alpha", "beta", "eta", "lambda", "deviance")

# Uncomment next line to save file
# save (items2export, file="objects_1MPE_Brazil.Rdata")



#####################
#### Mexico 2MPE #### 
#####################
rm(list=ls())

# Set working directory
savePathRep <- "/all_data/"

setwd (savePathRep)

all.processes  <- c("mexico_competing_final_result_sample_1a.Rda",  
                    "mexico_competing_final_result_sample_1b.Rda",
                    "mexico_competing_final_result_sample_2a.Rda",
                    "mexico_competing_final_result_sample_2b.Rda",
                    "mexico_competing_final_result_sample_3a.Rda",
                    "mexico_competing_final_result_sample_3b.Rda",
                    "mexico_competing_final_result_sample_4a.Rda",
                    "mexico_competing_final_result_sample_4b.Rda",
                    "mexico_competing_final_result_sample_5a.Rda",
                    "mexico_competing_final_result_sample_5b.Rda",
                    "mexico_competing_final_result_sample_6a.Rda",
                    "mexico_competing_final_result_sample_6b.Rda",
                    "mexico_competing_final_result_sample_7a.Rda",
                    "mexico_competing_final_result_sample_7b.Rda",
                    "mexico_competing_final_result_sample_8a.Rda",
                    "mexico_competing_final_result_sample_8b.Rda",
                    "mexico_competing_final_result_sample_9a.Rda",
                    "mexico_competing_final_result_sample_9b.Rda",
                    "mexico_competing_final_result_sample_10a.Rda",
                    "mexico_competing_final_result_sample_10b.Rda")


max.gelman.r <- c()
allChains <- mcmc.list()
for (k in 1:10){
  chain <- list()
  which.chains <- grep (paste0("sample_",k,"[ab]"), all.processes)
  
  # print sample iteration
  print (k)
  
  for (i in 1:length(which.chains)){
    load (all.processes[which.chains[i]])
    chain[[i]] <- as.mcmc(chainsComp); rm (chainsComp)
    # print process iteration within sample
    print (i)
  }
  
  ######################
  
  gelman.Comp <- gelman.diag (chain, multivariate=F)
  max.gelman <- max (gelman.Comp$psrf[,1], na.rm=T)
  print (max.gelman)

  allChains[[k]] <- as.mcmc.list (chain)
}

# Add convergence check for 20 chains all together
# This step ascertains that different chains *based on different samples*
# converge more or less on the same parameter values
Chains <- as.mcmc.list (list (
  allChains[[1]][[1]], allChains[[1]][[2]], 
  allChains[[2]][[1]], allChains[[2]][[2]], 
  allChains[[3]][[1]], allChains[[3]][[2]], 
  allChains[[4]][[1]], allChains[[4]][[2]], 
  allChains[[5]][[1]], allChains[[5]][[2]], 
  allChains[[6]][[1]], allChains[[6]][[2]], 
  allChains[[7]][[1]], allChains[[7]][[2]], 
  allChains[[8]][[1]], allChains[[8]][[2]], 
  allChains[[9]][[1]], allChains[[9]][[2]], 
  allChains[[10]][[1]], allChains[[10]][[2]]))

gelman.all <- gelman.diag (Chains, multivariate=F)
max.gelman <- max (gelman.all$psrf[,1], na.rm=T)
print (max.gelman)



theta <- delta <- beta <- alpha <- eta <- lambda <- c()
voteDis <- voteAgr <- voteUnc <- deviance <- c()
for (i in 1:10){
  process.1 <- allChains[[i]][[1]]
  process.2 <- allChains[[i]][[2]]
  
  th1 <- process.1[,grep("theta", colnames(process.1))]
  dt1 <- process.1[,grep("delta", colnames(process.1))]
  bt1 <- process.1[,grep("beta", colnames(process.1))]
  et1 <- process.1[,grep("^eta", colnames(process.1))]
  al1 <- process.1[,grep("^alpha", colnames(process.1))]
  lm1 <- process.1[,grep("lambda", colnames(process.1))]
  dis1 <- process.1[,grep("VoteDisagree", colnames (process.1))]
  agr1 <- process.1[,grep("VoteAgree", colnames (process.1))]
  unc1 <- process.1[,grep("probVoteUncertain", colnames (process.1))]
  dev1 <- process.1[,grep("deviance", colnames (process.1))]
  
  th2 <- process.2[,grep("theta", colnames(process.2))]
  dt2 <- process.2[,grep("delta", colnames(process.2))]
  bt2 <- process.2[,grep("beta", colnames(process.2))]
  et2 <- process.2[,grep("^eta", colnames(process.2))]
  al2 <- process.2[,grep("^alpha", colnames(process.2))]
  lm2 <- process.2[,grep("lambda", colnames(process.2))]
  dis2 <- process.2[,grep("VoteDisagree", colnames (process.2))]
  agr2 <- process.2[,grep("VoteAgree", colnames (process.2))]
  unc2 <- process.2[,grep("probVoteUncertain", colnames (process.2))]
  dev2 <- process.2[,grep("deviance", colnames (process.2))]
  
  th <- rbind (th1, th2)
  dt <- rbind (dt1, dt2)
  bt <- rbind (bt1, bt2)
  et <- rbind (et1, et2)
  al <- rbind (al1, al2)
  lm <- rbind (lm1, lm2)
  disagree  <- c(dis1, dis2)
  agree     <- c(agr1, agr2)
  uncertain <- c(unc1, unc2)
  dev  <- c(dev1, dev2)
  
  theta  <- rbind (theta, th)
  delta  <- rbind (delta, dt)
  beta   <- rbind (beta, bt)
  alpha  <- rbind (alpha, al)
  eta    <- rbind (eta, et)
  lambda <- rbind (lambda, lm)
  voteDis <- c(voteDis, disagree)
  voteAgr <- c(voteAgr, agree)
  voteUnc <- c(voteUnc, uncertain)
  deviance <- c(deviance, dev)
}

items2export <- list (theta, delta, alpha, beta, eta, lambda, voteDis, voteAgr, voteUnc, deviance)
names (items2export) <- c ("theta", "delta", "alpha", "beta", "eta", "lambda", "voteDis", "voteAgr", "voteUnc", "deviance")

# Uncomment next line to save file
# save (items2export, file="objects_2MPE_Mexico.Rdata")


#####################
#### Mexico 1MPE #### 
#####################
rm(list=ls())

# Set working directory
savePathRep <- "/all_data/"
setwd (savePathRep)

all.processes <- c("mexico_pressure_final_result_sample_1a.Rda",
                   "mexico_pressure_final_result_sample_1b.Rda",
                   "mexico_pressure_final_result_sample_2a.Rda",
                   "mexico_pressure_final_result_sample_2b.Rda",
                   "mexico_pressure_final_result_sample_3a.Rda",
                   "mexico_pressure_final_result_sample_3b.Rda",
                   "mexico_pressure_final_result_sample_4a.Rda",
                   "mexico_pressure_final_result_sample_4b.Rda",
                   "mexico_pressure_final_result_sample_5a.Rda",
                   "mexico_pressure_final_result_sample_5b.Rda",
                   "mexico_pressure_final_result_sample_6a.Rda",
                   "mexico_pressure_final_result_sample_6b.Rda",
                   "mexico_pressure_final_result_sample_7a.Rda",
                   "mexico_pressure_final_result_sample_7b.Rda",
                   "mexico_pressure_final_result_sample_8a.Rda",
                   "mexico_pressure_final_result_sample_8b.Rda",
                   "mexico_pressure_final_result_sample_9a.Rda",
                   "mexico_pressure_final_result_sample_9b.Rda",
                   "mexico_pressure_final_result_sample_10a.Rda",
                   "mexico_pressure_final_result_sample_10b.Rda")
                       
                       

# Import converged processes to check graphs
max.gelman.r <- c()
allChains <- mcmc.list()
for (k in 1:10){
  chain <- list()
  which.chains <- grep (paste0("sample_",k,"[ab]"), all.processes)
  
  # print sample iteration
  print (k)
  
  for (i in 1:length(which.chains)){
    load (all.processes[which.chains[i]])
    chain[[i]] <- as.mcmc(chainsPress); rm (chainsPress)
    # print process iteration within sample
    print (i)
  }
  
  gelman.Comp <- gelman.diag (chain, multivariate=F)
  max.gelman <- max (gelman.Comp$psrf[,1], na.rm=T)
  print (max.gelman)
  max.gelman.r <- c(max.gelman.r, max.gelman)
  
  allChains[[k]] <- as.mcmc.list (chain)
}

# Add convergence check for 20 chains all together
# This step ascertains that different chains *based on different samples*
# converge more or less on the same parameter values
Chains <- as.mcmc.list (list (
                        allChains[[1]][[1]], allChains[[1]][[2]], 
                        allChains[[2]][[1]], allChains[[2]][[2]], 
                        allChains[[3]][[1]], allChains[[3]][[2]], 
                        allChains[[4]][[1]], allChains[[4]][[2]], 
                        allChains[[5]][[1]], allChains[[5]][[2]], 
                        allChains[[6]][[1]], allChains[[6]][[2]], 
                        allChains[[7]][[1]], allChains[[7]][[2]], 
                        allChains[[8]][[1]], allChains[[8]][[2]], 
                        allChains[[9]][[1]], allChains[[9]][[2]], 
                        allChains[[10]][[1]], allChains[[10]][[2]]))

gelman.all <- gelman.diag (Chains, multivariate=F)
max.gelman <- max (gelman.all$psrf[,1], na.rm=T)
print (max.gelman)



theta <- delta <- beta <- alpha <- eta <- lambda <- deviance <- c()
for (i in 1:10){
  process.1 <- allChains[[i]][[1]]
  process.2 <- allChains[[i]][[2]]
  
  th1 <- process.1[,grep("theta", colnames(process.1))]
  dt1 <- process.1[,grep("delta", colnames(process.1))]
  bt1 <- process.1[,grep("beta", colnames(process.1))]
  et1 <- process.1[,grep("^eta", colnames(process.1))]
  al1 <- process.1[,grep("^alpha", colnames(process.1))]
  lm1 <- process.1[,grep("lambda", colnames(process.1))]
  dev1 <- process.1[,grep("deviance", colnames (process.1))]
  
  th2 <- process.2[,grep("theta", colnames(process.2))]
  dt2 <- process.2[,grep("delta", colnames(process.2))]
  bt2 <- process.2[,grep("beta", colnames(process.2))]
  et2 <- process.2[,grep("^eta", colnames(process.2))]
  al2 <- process.2[,grep("^alpha", colnames(process.2))]
  lm2 <- process.2[,grep("lambda", colnames(process.2))]
  dev2 <- process.2[,grep("deviance", colnames (process.2))]
  
  th <- rbind (th1, th2)
  dt <- rbind (dt1, dt2)
  bt <- rbind (bt1, bt2)
  et <- rbind (et1, et2)
  al <- rbind (al1, al2)
  lm <- rbind (lm1, lm2)
  dev  <- c(dev1, dev2)
  
  theta  <- rbind (theta, th)
  delta  <- rbind (delta, dt)
  beta   <- rbind (beta, bt)
  alpha  <- rbind (alpha, al)
  eta    <- rbind (eta, et)
  lambda <- rbind (lambda, lm)
  deviance <- c(deviance, dev)
}

items2export <- list (theta, delta, alpha, beta, eta, lambda, deviance)
names (items2export) <- c ("theta", "delta", "alpha", "beta", "eta", "lambda", "deviance")

# Uncomment next line to save file
# save (items2export, file="objects_1MPE_Mexico.Rdata")



#####################################################################
#### Checking convergence for chains based on mean sample values ####
#####################################################################

#####################
#### Brazil 2MPE #### 
#####################
rm(list=ls())

# Set working directory
savePathRep <- "/all_data/"
setwd (savePathRep)

# Change names in the next two lines
all.processes <- c("brazil_competing_final_result_sample_meana.Rda",
                   "brazil_competing_final_result_sample_meanb.Rda")


max.gelman.r <- c()
allChains <- mcmc.list()
chain <- list()
which.chains <- grep (paste0("sample_mean","[ab]"), all.processes)
  
for (i in 1:length(which.chains)){
  load (all.processes[which.chains[i]])
  chain[[i]] <- as.mcmc(chainsComp); rm (chainsComp)
    # print process iteration within sample
}
  
print (max (gelman.diag (chain, multivariate=F)$psrf[,1], na.rm=T))


#####################
#### Brazil 1MPE #### 
#####################
rm(list=ls())

# Set working directory
savePathRep <- "/all_data/"
setwd (savePathRep)

# Change names in the next two lines
all.processes <- c("brazil_pressure_final_result_mean_a.Rda",
                   "brazil_pressure_final_result_mean_b.Rda")


max.gelman.r <- c()
allChains <- mcmc.list()
chain <- list()
which.chains <- grep (paste0("mean_","[ab]"), all.processes)

for (i in 1:length(which.chains)){
  load (all.processes[which.chains[i]])
  chain[[i]] <- as.mcmc(chainsPress); rm (chainsPress)
  # print process iteration within sample
}

print (max (gelman.diag (chain, multivariate=F)$psrf[,1], na.rm=T))


#####################
#### Mexico 2MPE #### 
#####################
rm(list=ls())

# Set working directory
savePathRep <- "/all_data/"
setwd (savePathRep)

# Change names in the next two lines
all.processes <- c("mexico_competing_final_result_sample_meana.Rda",
                   "mexico_competing_final_result_sample_meanb.Rda")


max.gelman.r <- c()
allChains <- mcmc.list()
chain <- list()
which.chains <- grep (paste0("sample_mean","[ab]"), all.processes)

for (i in 1:length(which.chains)){
  load (all.processes[which.chains[i]])
  chain[[i]] <- as.mcmc(chainsComp); rm (chainsComp)
  # print process iteration within sample
}

print (max (gelman.diag (chain, multivariate=F)$psrf[,1], na.rm=T))


#####################
#### Mexico 1MPE #### 
#####################
rm(list=ls())

# Set working directory
savePathRep <- "/all_data/"
setwd (savePathRep)

# Change names in the next two lines
all.processes <- c("mexico_pressure_final_result_sample_meana.Rda",
                   "mexico_pressure_final_result_sample_meanb.Rda")


max.gelman.r <- c()
allChains <- mcmc.list()
chain <- list()
which.chains <- grep (paste0("mean","[ab]"), all.processes)

for (i in 1:length(which.chains)){
  load (all.processes[which.chains[i]])
  chain[[i]] <- as.mcmc(chainsPress); rm (chainsPress)
  # print process iteration within sample
}

print (max (gelman.diag (chain, multivariate=F)$psrf[,1], na.rm=T))


